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Anisotropy of the permeability tensor in statistically uniform porous media of sizes used in typical 
computer simulations is studied. Although such systems are assumed to be isotropic by default, we 
show that de facto their anisotropic permeability can give rise to significant changes of transport 
parameters such as permeability and tortuosity. The main parameter controlling the anisotropy 
is a/L, being the ratio of the obstacle to system size. Distribution of the angle a between the 
external force and the volumetric fluid stream is found to be approximately normal, and the standard 
deviation of a is found to decay with the system size as (a/L) d ^ 2 , where d is the space dimensionality. 
These properties can be used to estimate both anisotropy-related statistical errors in large-scale 
simulations and the size of the representative elementary volume. 

PACS numbers: 47.56.+r,47.15.G-,91.60.Np 



I. INTRODUCTION 

A standard method of modeling a uniform, isotropic 
porous medium (e.g. a column of sand) is to place ran- 
domly many identical objects that are impermeable to 
fluid (e.g. solid spheres) in an initially empty volume [1- 
8]. Since the objects are placed uniformly in the whole 
system, one might expect that randomness in their exact 
locations is irrelevant in the sense that the bulk volu- 
metric fluid stream will be parallel to the external force 
(e.g. gravitation). This would be the case if the sys- 
tem was large enough. However, in computer simula- 
tions and in artificial laboratory systems (used in par- 
ticle image velocimetry measurements [9]), usually rela- 
tively small systems are utilized that contain at most a 
few thousands of "grains" — far less than billions of sand 
grains in a typical experimental setup. Since randomly 
distributed grains tend to form channels of random orien- 
tations, small porous systems are very sensitive to local 
fluctuations of the grain distribution. Under such con- 
ditions the direction of the volumetric fluid stream can 
differ significantly from that of the external force. Con- 
sequently, a system that was supposed to be isotropic, 
may de facto be rather highly anisotropic. The aim of 
this paper is a detailed analysis of this phenomenon in a 
two-dimensional (2D) flow. 

A porous medium is anisotropic to flow if the perme- 
ability tensor is anisotropic. Usually anisotropy of the 
permeability tensor is either assumed explicitly [10] or it 
appears naturally as an expected consequence of a micro- 
scopic model [11-14]. In the former case one works en- 
tirely on a macroscopic level, whereas the latter approach 
tries to connect the observed macroscopic anisotropy of 
real porous materials with their microscopic geometry 
and structure. Permeability anisotropy caused by a fi- 
nite size of a model system has not been regarded as 
an important factor so far, although some research tech- 



niques, e.g. numerical simulations, concentrate on artifi- 
cially small porous systems. The reason for this lies in the 
fact that numerical flow simulations in complex porous 
geometries are extremely tedious and require extensive 
computer power and resources. A common strategy has 
been to perform calculations for just a few systems that 
are as large as possible [7, 15] . In contrast to this, here we 
solve the flow equations for hundreds or even thousands 
of different porous systems of small to medium sizes and 
then extrapolate the results to the limit of an infinitely 
large system. This method was already used in [6] to 
detect a small, systematic deviation of the flow tortu- 
osity from several theoretical formulas, with an ad hoc 
interpretation of this phenomenon as a consequence of 
the finite-size anisotropy. Therefore, in this paper we 
present a systematic study of finite-size anisotropy in 
a two-dimensional model of statistically uniform porous 
media. 

The structure of the paper is as follows. Section II 
specifies the model and the numerical techniques used. 
Main results are provided in Sec. III. Next, in Section IV 
we develop a simple theory to account for the asymptotic 
behavior of the angle between the external force and the 
volumetric fluid flux. Finally, the results are discussed in 
Sec. V. 



II. MODEL 

In this study we use a model of freely overlapping 
squares [2, 6, 16]. In essence, this is a two-dimensional 
lattice model with a porous matrix modelled as a union 
of freely overlapping identical solid squares of size a x a 
lattice units (l.u.) placed uniformly at random locations 
in a square lattice L x L l.u. (1 < o < L). The squares 
are fixed in space but free to overlap, and their sides co- 
incide with the underlying lattice. The remaining void 
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space is filled with a fluid. A constant, external force is 
imposed on the fluid to model the gravity and we allow 
an angle (/3) between the force and the system side to 
be arbitrary (note that in [2, 6, 16] only the case (3 = 
was considered). Periodic boundary conditions are im- 
posed in both directions to minimize finite-size effects. 
The porosity (<ft) is calculated as the ratio of unoccupied 
lattice nodes to the entire system volume (L 2 ). The flow 
equations are solved in the creeping flow regime using 
the Lattice Boltzmann Model (LBM) [17] with a single 
relaxation time collision operator [18] (see [6] for imple- 
mentation details). 

The model has three adjustable parameters: </>, a, and 
L. The first one corresponds directly to the macroscopic 
porosity. The value of a affects the percolation threshold 
4> c , which is a decreasing function of a from 4> c w 0.4073 
(the standard site percolation threshold, a = 1) [19] to 
4> c ~ 0.3333 (the continuous percolation threshold of 
aligned squares, a — > oo) [20]. As the model is solved 
using the LBM method without a numerical grid refine- 
ment [6], the minimum value of a is 4 (this is the min- 
imum length scale for the LBM method to resolve the 
macroscopic Navier-Stokes equations [17]). The value of 
L controls the finite-size effects through the dimension- 
less ratio a/L, which should be as small as possible to 
mimic an infinite system. 

Anisotropy of fluid flow in the above-defined model will 
be investigated through Darcy's law [10] 



q = Kg, 



(1) 



where q is the volumetric fluid flux, K is a symmetric 
tensor of the hydraulic conductivity, and g is the unit 
vector in the direction of the gravitational field. 



III. NUMERICAL RESULTS 

A. Tests on K 

The basic concept of transition from microscopic laws 
of hydrodynamics to macroscopic laws of transport in 
porous media is the representative elementary volume 
(REV), i.e. the smallest volume such that a measurement 
over it will yield a value representative of the whole [10]. 
Darcy's law (1) is, in principle, applicable only to sys- 
tems that are larger than an REV, whereas significant 
anisotropy is expected in systems smaller than an REV. 
Hence, the primary question is whether or not Eq. (1) 
can be used to study anisotropy in small-size systems. 
To answer this we performed several simulations on K , 
with its elements computed from q(g) for g = x and 

g = y- 

First, the symmetry of K was examined by quantifying 
the value of a dimcnsionless parameter given by 
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FIG. 1: Streamlines through the same porous system (L = 
100 l.u., a — 4 l.u., 4> = 0.7) for two different alignments of 
the external force g. The grey squares represent the solid 
part of the medium, and the remaining space is open to fluid 
flow. Left panel: g is parallel to the x axis (f3 = 0) and the 
specific discharge q makes an angle a m 21° with x. Right 
panel: (5 w —22° (calculated from Eq. (3)) and the angle 
between the specific discharge q and the x axis is a « 0.7°. 
For the ease of display, two selected streamlines and their 
counterparts in both panels are given in color. 



Furthermore, by choosing L = 100 l.u., a = 4 l.u., and 
4> = 0.7, 0.9, eighty different (i.e., fourty systems for each 
4>) statistically uniform porous systems were constructed, 
for which e < 0.5% was found. This ensures that K is 
symmetric within 0.5% numerical errors in its elements. 
In a subsequent analysis we enforced _K"-symmetry via re- 
placing its off-diagonal elements (K xy and K yx ) by their 
arithmetic mean, which ensures that K is diagonizable. 

Second, the tensorial properties of K were examined 
by checking whether Eq. (1) can be used for an arbitrary 
g. In particular, this equation predicts that if the mean 
flow direction (q) is aligned with the x-axis, the angle 
(3 between the external force (g) and the n-axis should 
satisfy 



K 

tan^ = 



K 



(3) 



.'/// 



This relation was examined for several systems, of which 
one is shown in Fig. 1, wherein, two streamline patterns 
for the same system (L = 100 l.u., a — 4 l.u., <fi = 0.7) 
with different g-orientations are visualized. In the left 
panel, the external force is parallel to the (horizontal) 
x-axis (g = x), resulting in an angle of a w 21° between 
the vector of the specific discharge (q) and the x axis. 

In the right panel, a force of the same magnitude makes 
an angle f3 w —22° computed from (3); as expected, the 
angle between q and the x axis practically vanishes (a ~ 
-0.7°). 



B. 



Tests on a 



(2) 



A natural measure of anisotropy for a particular porous 
system is the angle between the vectors g and q. As this 
angle depends on the orientation of g, following standard 
procedures in computer simulations, we fix g = x. We 
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FIG. 2: <r Q = yj (a 2 ) as a function of L for a — 4 l.u. and 
two porosities <j> = 0.7, 0.9, with error bars at 95% confidence 
level. The lines represent fits to the power law a a oc L _1 for 
L > f00 l.u. 

verified that in this case the numerical value of a (angle 
between q and g) satisfies (a) w 0, which follows from 
symmetry arguments, and then calculated 

o a = y^y. (4) 

In the above equation, (• • • ) denotes an average over dif- 
ferent random porous systems. The results for a = 4 l.u., 
4> = 0.7, 0.9, and several system lengths L are shown in 
Fig. 2. For L > 100 l.u. the data were fitted to 

a a oc L~ 5 , (5) 

which yielded 5 w 0.96(6) for = 0.7 and 5 w 1.00(3) for 
<f> = 0.9. This suggests 5=1, i.e. 

CTaOcL -1 , L>1. (6) 

This relation does not hold for small L (L < 50 in Fig. 2), 
for some realizations of such systems are likely to exhibit 
extreme anisotropy with a so large that sin a cannot be 
approximated by a (for L = 50 the angle between q and 
g can be as large as 45°). 

Next we investigated statistical distribution of a- values 
in different random systems with fixed L, a, and <f>. In all 
cases this distribution closely resembled the normal dis- 
tribution N(0,a^). Qualitative verification of this con- 
jecture is presented in Fig. 3, which depicts the empirical 
cumulative distribution function (CDF) for two different 
system sizes L (small symbols) together with the corre- 
sponding theoretical CDFs of the normal distribution, 

F(a) = 1 + CTf( ^ ) . (7) 

The numerical data are in good agreement with (7). 
A quantitative comparison of the a-distribution with 
iV(0,<7^) was performed using the Kolmogorov-Smirnov 
test (at confidence level 95%). Out of all data points 
shown in Fig. 2 only the one corresponding to L = 50 
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FIG. 3: Cumulative distribution function (CDF) of a for (f> — 

0. 7, a = 4 l.u. and two system sizes L — 100 and L — 400 

1. u. (dots). The empirical CDF was determined using 500 (for 
L — 100) and 80 (L = 400) numerical samples. Solid lines 
represent theoretical CDF of the normal distribution, Eq. (7), 
with G a = \/ {oc 2 )- 



and <fi = 0.7 did not pass the test, which in part is due to 
extremely large number of different samples (2000) used. 

As mentioned before, the value of a determines the 
percolation threshold <j> c for small porosities, and can be 
considered as a relevant parameter independent of (j> and 
L. For porosities much larger than <f> c , however, the con- 
nectedness and overlapping of individual randomly gener- 
ated solid squares becomes irrelevant. In this case, using 
scaling arguments, one can expect that <fi and a/L are 
the only relevant parameters. Mathematically, this can 
be formulated as a similarity ansatz: 

a a (a,L,4>) K,ty(a/L,4>), r>»0 c ,i>a. (8) 

where ^ is a similarity function. A direct numerical ver- 
ification of this conjecture is difficult, as it requires av- 
eraging over many independent samples, which is rather 
a time-consuming job for large a. Instead of this, we 
concentrated on a single parameter set with <f> — 0.7 and 
a/L = 0.04 that led to results demonstrated in Fig. 4 
shown as cross symbols. These data were fitted to an 
ad-hoc formula 

a a (a) = ci + c 2 cxp(-a/c 3 ) (9) 

with three adjustable parameters c\, C2, and c 3 . The 
best-fit value of c 3 w 0.6 indicates that the approximation 
(8) can be safely used for a > 4. 

Finally, we investigated the dependence of o~ a on poros- 
ity. One expects that a a should decrease from ss 45° at 
the percolation threshold <j> c (a single, randomly oriented 
conducting channel) to 0° at <p — 1 (completely perme- 
able system). As shown in Fig. 5, our numerical results 
generally agreed with this picture. However, a a did not 
converge to its limiting value as 4> — > 1. Instaed, it 
saturated at a positive value, which is independent of 4>. 
Due to this rather unexpected result, we ensured that 



FIG. 4: a a as a function of a for cj> = 0.7 and a/L — 0.04 
(x symbols). The data come from N = 500 independent 
porous systems for a < 4 l.u. and N = 200 for a > 4 l.u. The 
error bars were calculated at the 95% confidence level. The 
solid line represents the best fit to Eq. (9) . 
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FIG. 5: (J a (in degrees) as a function of porosity (j> for L — 100 
l.u. Inset: the product a a L (deg. x l.u.) for L = 100 l.u. 
(cross symbols) and L — 200 l.u. (circle symbols). All data 
obtained for a = 4 l.u.; error bars at 95% confidence level 
obtained from 200 independent porous systems. The arrows 
show the percolation threshold 4> c ~ 0.417. Filled circles show 
theoretical values of a a for (f> — 4> c , L — > oo (45°) and 4> — 1 
(0°). 



neither discretization errors nor large relaxation times 
affect the numerical data obtained for large porosities. 
We also verified that the system size used, L = 100, is 
sufficiently large for relation (6) to hold. This is clearly 
seen in the inset of Fig. 5, which depicts the product a a L 
for L = 100 and L = 200. The data for different L col- 
lapsed in a broad range of > 0.55. Porosities less than 
« 0.55 are in a vicinity of the percolation critical point, 
at which <r a is expected to converge to 45° as L — > oo, 
and hence the product a a L should diverge at 4> c . As the 
system size L — > oo, it is possible that the the porosity 
range, for which scaling relation (6) does not hold, di- 
minishes according to a power law. This behavior is a 
typical finite-size effect near a critical point [21]. 

To explore the reason why <r a docs not tend smoothly 
to as approaches 1, we inspected the streamlines in 
high-porosity systems exhibiting large anisotropy. An ex- 



FIG. 6: Streamlines in a high-porosity system (4> = 0.95) 
with L = 100 l.u. and a = 4 l.u. for two different flow types: 
(a) hydrodynamic (a £s 15°); (b) electric (a £s 0.6°). The 
electric flow was calculated from a solution of the Laplace 
equation with periodic boundary conditions and a constant 
electric field parallel to the a;-axis. Note that the distribution 
of obstacles and the orientation of an external body force is 
identical in both panels. 



treme example of such a system, generated with cf) — 0.95, 
is shown in Fig. 6a. At this high porosity, overlapping of 
individual obstacles is negligible, and the solid part of the 
system is made up of separate islands (that could corre- 
spond, for example, to a cross-section of a porous medium 
made of parallel fibers [22, 23]). Because the obstacles 
were placed uniformly and randomly in the whole system, 
their local concentration varies, and they form several 
larger groups of obstacles with relatively small distances 
between group members. Since fluid flux through a 2D 
channel is proportional to its width squared, most of the 
fluid flow takes place in relatively wide 'channels' between 
the groups. In other words, owing to no-slip boundary 
conditions on the obstacle surfaces, the fluid passes most 
easily in the inter-connected regions of low local obstacle 
concentration (high local porosity), whereas the regions 
of high local obstacle concentration (low local porosity) — 
even if occupied by separate obstacles — act effectively as 
large, almost impenetrable barriers. This solid-fluid 're- 
pulsion' effect is not present in electric current flows (for 
the current intensity is proportional to the first power 
of a conductive channel width). For this reason, a high- 
porosity system which is highly anisotropic to fluid flow 
(a ~ 15° in Fig. 6a) exhibits a marginal anisotropy to 
electric current flow (|a| < 1°) as depicted in Fig. 6b. 



C. Tests on principal values 

Mathematically, a porous system is anisotropic to flow 
if and only if at least two of the principal values of K are 
diffrcnt. In the present case K has two eigenvalues (prin- 
cipal permeabilities) K + and K ~ which can be ordered 
such that K + > K~ . Their ratio, 
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FIG. 7: Cumulative distribution function of r = K~ / K + for 
(j> = 0.7, a = 4 l.u., L = 100 l.u., calculated using AT = 340 
different porous systems. Dashed line represents the best fit 
to the CDF of the normal distribution with (r) w 0.75 and 
o r « 0.12. 

is equivalent to the ratio of the minimum to maximum 
permeabilities of a given porous system, and hence is a 
proper measure of its anisotropy [11]. The more r devi- 
ates from 1, the more anisotropic the system is. 

We first verified that, as expected, the angle between 
the main principal axis and the x-axis was distributed 
uniformly in the range (— 7r/2,7r/2] (data not shown). 
Then the CDF of r was determined for a particular case 
with (j) — 0.7, a = 4 l.u., and L = 100 l.u. As can be 
seen in Figure 7, the distribution of r can be quite well 
fitted to the normal distribution iV(0.75, 0.12 2 ). How- 
ever, this is only an approximation, as in the present 
case CDF(r) = 1 for r > 1. 

IV. DISTRIBUTION OF a FOR LARGE L 

Consider a porous system of size L x L l.u. subject 
to an external force along the x axis. Let Ar denote 
the total displacement of a fluid particle as it passes the 
system between the opposite boundaries. While the x- 
componcnt of Ar is a constant (equal to the system size 
L), the y component (which we shall call lateral displace- 
ment and denote Ay) varies for different streamlines. If 
we calculate the average (Ay) over all fluid particles, the 
angle a between the volumetric fluid flux q and the x 
axis will satisfy 

tan a =<M. (11) 

hi 

If a is sufficiently small, this equation simplifies to 




FIG. 8: A two-dimensional porous system (a), which can be 
regarded as a two-layer system perpendicular to the external 
force (b) or a group of four smaller subsystems (c) . 

as two subsystems of size 2L x L l.u. (labelled A, B) 
or four subsystems of size L x L l.u. (labelled 1,2,3,4). 
Each of the small subsystems has its own permeability 
tensor Kj, volumetric fluid flux qj, angle ctj between 
the x axis and qj, and mean lateral displacement (Ayj) 
with j — 1, ... ,4. Since the distribution of obstacles is 
uniform, porosities of each small subsystem is approx- 
imately equal to </>, and the mean lateral displacements 
(Ayj) can be considered as independent random variables 
drawn from the same distribution. Subsystems 1 and 2 
form layer A orthogonal to the external force. One may 
expect that the fluid streams passing through subsystems 
1 and 2 are approximately the same in magnitude, and 
so the mean lateral displacement of the fluid, as it passes 
through layer A, can be approximated by 

( ^ m <M+£b> p,, 

Similarly, the mean lateral displacement in the layer 
B can be approximated by (Ays) ~ ((Ay$) + (Aj/4))/2. 
The mean lateral displacements of the fluid in layers A 
and B are practically independent of each other. This can 
be justified by an example of soil made of several hori- 
zontal and anisotropic layers — in this case the mean flow 
direction in a layer will depend only on the permeability 
tensor of this layer. This implies that the total lateral 
displacement of the fluid in the whole system ({Ay)) is 
approximately given by 

(Ay) « (Ay A ) + (Ay B ) « \ ]T(A % ). (14) 

i=i 

If L is large, then a becomes sufficiently small for ap- 
proximation (12) to be valid. In this case Eqs. (14) and 
(12) lead to 

1 4 

where a is calculated for the whole, 2L x 2L system. As- 
suming that a, are independent random variables drawn 
from the same distribution with mean 0, one arrives at 



Let us consider a porous system of size 2L x 2L l.u. and 
porosity <p. As shown in Fig. 8, it can also be regarded 



<r a (2L) w -a a (L), 



(16) 



6 



which immediately leads to (6). 

Equation (15) can be used iteratively to obtain 

a(2 fc L)«7fc5>*( L )' fc = l,2,... 



10 



(17) 



where the arguments of a and ctj (i.e. 2 fc L and L) in- 
dicate the system size. The right-hand side of this for- 
mula is an arithmetic mean of independent random vari- 
ables with finite mean and variance, and — due to the 
central limit theorem — converges to normal distribution 
as k — > oo. This explains why the distribution of a for a 
sufficiently large system size L can be approximated by 
a normal distribution (see Fig. 3). 

The above can be readily extended to flows in an arbi- 
trary space dimension d. We skip the details and report 
only the final conclusions. First, 



ozL- s , 5 = d/2 



(18) 



of 



for sufficiently large L. Second, the distribution 
tends to the normal distribution as L — > oo. 

Equation (18) implies that anisotropy effects diminish 
with system size most quickly in three-dimensional (3D) 
systems (a a oc L~ 3 / 2 ). Note, however, that the most im- 
portant factor in computer simulations is the total num- 
ber of lattice nodes (or volume) V in the system. Using 
this quantity, equation (18) can be written as 



oc V 



-1/2 



(19) 



irrespective of d. Thus, anisotropy of the permeability 
tensor should be equally important (and difficult to ac- 
count for) in computer simulations carried out in any 
space dimension. 

It is important to verify Eq. (18) for space dimensions 
d ^ 2. While at the moment our software cannot be used 
for d = 3, the case d = 1 can be tackled by studying a 
quasi one-dimensional system of size K x L with fixed 
K and L — ► oo. Note that in this case Eq. (18) should 
hold irrespective of whether the longer side of the system 
is parallel or perpendicular to the external force. The 
results, obtained for a = 4 l.u., (f) = 0.7, K = 100, and 
L ranging from 100 to 800 l.u. are shown in Fig. 9 and 
confirm the validity of Eq. (18). 

Equations (8) and (18) allow to factorize <r a (a, L, <j>): 



(20) 



where $ is a function. This relation can be expected to 
hold in general only if L 3> a and <f> is sufficiently far away 
from the critical porosity <p c . In a general case, a is to be 
interpreted as a characteristic system length (such as the 
diameter of discs, in case the porous matrix is made of 
discs rather than squares), and $ depends on the system 
in question. 
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FIG. 9: a a for a quasi one-dimensional system of size K x 
L with K fixed at 100 l.u. and L growing from 100 to 800 
l.u. for channel axis parallel (x) and perpendicular (o) to 
the external force (a = 4 l.u., <j> = 0.7, error bars at 95% 
confidence level). Dashed line represents a single fit to a a oc 
1/y/L for all data points. 



V. DISCUSSION AND CONCLUSIONS 

Our results show that permeability anisotropy in sta- 
tistically uniform porous systems of sizes typically used 
in computer simulations is a significant factor. The main 
parameter controlling this phenomena, especially at high 
porosities, is the ratio a/L. For the model considered 
here, the asymptotic regime is observed for a/L < 0.04. 
In this regime the distribution of the angle a between the 
external force and the volumetric fluid flux is very close 
to Gaussian, with the standard deviation diminishing as 
{a/L) d ' 2 . 

Although this conclusion is based on numerical results 
obtained for a particular model of a two-dimensional flow, 
it is expected to apply to a wide class of porous systems 
with randomly distributed identical solid matrices, such 
as squares, disks or spheres. This observation can be 
used to estimate the anisotropy-related statistical error 
in large-scale simulations, where often only one large sys- 
tem is considered for each parameter set [15]. To this 
end it is enough to perform many independent simu- 
lations in small- and medium-size systems, verify that 
a a oc (a/L) d / 2 , and extrapolate a a (L) to the required 
value of L. Next, assuming that the distribution of a is 
normal, one obtains the complete information about the 
error related to the anisotropy of the permeability tensor. 

Magnitude of permeability anisotropy could serve as 
a good indicator of how the size of a model system 
compares with that of a REV. We found that even for 
a/L = 0.04 the angle between the external force and 
the volumetric fluid flux can be as large as 20°, and the 
permeability can vary with the orientation of the ex- 
ternal force by a factor of 2. The value below which 
the anisotropy effects are small enough to be practically 
negligible is a/L rts 0.01, as in this case a a < 2°, i.e. 
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| or j < 6° with probability p k, 0.99. This enables to esti- 
mate the size of a REV in the model considered here as 
w 400 x 400 l.u. 

It is interesting to note that most of the simulations 
carried out so far for 2D systems do not meet the criterion 
of a/L < 0.01, mainly because they used models with 
large a. In previous studies on two-dimensional flows 
in various statistically uniform porous media, many re- 
searchers used a/L- values ranging from 0.02 [24], through 
0.026 [4], 0.04 [1, 5], 0.05 [6, 25] to 0.1 [2, 3, 16], usually 
assuming their systems to be isotropic. In view of our 
present findings, validity of this assumption in some of 
these cases is questionable and requires verification. Gen- 
erally, one should expect that the threshold value of a/L 
below which the permeability anisotropy is negligible is 
not universal, but depends on the geometry and struc- 
ture of the porous medium, especially on its porosity and 
space dimensionality. 

Anisotropy is a phenomenon independent of the 
boundary conditions. Periodic boundary conditions used 
in this paper facilitate measurement of the permeability 
tensor and reduce finite-size (boundary) effects. Other 
boundary conditions could mask, but would not elimi- 
nate anisotropy effects. For example, using solid walls 
along the fluid flow would fix the orientation of the 
fluid stream, however, the system would respond to such 
boundaries with an internal pressure gradient [10], which 



would change (and complicate measurement of) the ori- 
entation of the effective force acting on the fluid. 

Finite-size permeability anisotropy in three- 
dimensional small porous systems remains an open 
problem. Typical system sizes used in numerical 3D 
simulations are L w 100 l.u. The ratio a/L is thus much 
larger in 3D than in 2D simulations and ranges from 
0.06 [8, 26], through 0.1 [25], 0.125 [15, 27], to 0.33 [4]. 
The magnitude of permeability anisotropy is usually 
neglected. One exception is the paper by Verberg and 
Ladd [4] , who published the off-diagonal elements of the 
permeability tensor. Their data for a single configuration 
of randomly distributed spheres suggests that a a is a 
decreasing function of the porosity and varies from 
<j a w 3° for <f> = 0.464 to a a w 18° for = 0.087. This is 
in agreement with our present findings for a 2D system 
and indicates that permeability anisotropy is especially 
important close to the percolation threshold. 
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